% Evaluate total compliance by Pulse Pressure Method 

aantal=length(pressure);

meanP=mean(pressure);
meanQ=mean(flow);

TVR=meanP/meanQ;
CO=meanQ*60/1000;
T=max(time);
HR=60/T;


% compute input impedance

Zin=fft(pressure)./fft(flow);
Zinmod=abs(Zin);
Zinphase=unwrap(angle(Zin));
freq=(1/T)*((1:aantal)-1);


% THE PULSE PRESSURE METHOD


C_init=1.5;						% initial compliance value in ml / mmHg
error=0.001;						% max error 0.001%


% true PP is defined as the max - min pressure in early systole

[bps ibps]=max(pressure(1:aantal));
%truePP=max(pressure(1:aantal))-min(pressure(1:ibps));			% mmHg	
truePP=max(pressure(1:aantal))-min(pressure(1:aantal));			% mmHg	

%substract DC level and do FFT on flow


Qtofft=flow-meanQ;

fftQ=fft(Qtofft);
magQ=2*abs(fftQ)/length(fftQ);
phaseQ=unwrap(angle(fftQ));


C_WK2=C_init;
fftP_WK2(1)=0;
comperr=2*error;


% application PPM: compute P-response of 2-el WK model to flow

while abs(comperr)>error
	 Z_WK2(1)=TVR;
	 for j=2:aantal
	    freq(j)=(j-1)/T;
	    Z_WK2(j)=TVR/(1+i*2*pi*freq(j)*TVR*C_WK2);
	    fftP_WK2(j)=Z_WK2(j)*fftQ(j);
	 end

	P_WK2=2*real(ifft(fftP_WK2(1:aantal/2),aantal))+meanP;

 	PP_WK2=max(P_WK2)-min(P_WK2);
 	comperr=100*(truePP-PP_WK2)/truePP;
	C_WK2=C_WK2-comperr*C_WK2/100;
end;

fig1=figure;

subplot(2,1,1);
title('input impedance');
xlabel('frequency (Hz)');
ylabel('impedance (mmHg/(ml/s))');
semilogy(freq(1:11),abs(Zin(1:11)),'m+',freq(1:11),abs(Z_WK2(1:11)),'y');
legend('measured','model');

subplot(2,1,2);
title('pressure');
xlabel('time (s)');
ylabel('pressure (mmHg)');
plot(time,pressure,time,P_WK2);
legend('measured','model');
fprintf(1,'estimated C (PPM): %g\n',C_WK2)

